\name{smith}
\alias{smith}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{Smith Form (approximate algorithm)}
\description{
%%  ~~ A concise (1-5 lines) description of what the function does. ~~
}
\usage{
smith(A, tol)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
  \item{A}{Polynomial matrix.}
  \item{tol}{Tolerance for the approximate algorithm (real).}
}
\details{
This function applies the approximate algorithm in [1] to matrix A. Thus, if A is a consistent
 estimate of a certain matrix A0, the output of the function is a consistent estimate of A0.
}
\value{
%%  ~Describe the value returned
%%  If it is a LIST, use
%%  \item{comp1 }{Description of 'comp1'}
%%  \item{comp2 }{Description of 'comp2'}
%% ...
}
\references{
[1] Arbues, I., Ledo, R. Estimation of the Smith form and generalized VEC models.
(to appear).}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
%% ~~objects to See Also as \code{\link{help}}, ~~~
}
\examples{
x=polynomial(coef=c(0,1))
A=PolyMatrix(c(2,2),x^2,1+x^3,1-x^2,x^2)
smith(A,0.1)

## The function is currently defined as
function (A, tol) 
{
    rho = 1 + tol
    tol2 = 1e-10
    n = dim(A)[1]
    m = dim(A)[2]
    U = unitPM(m)
    V = unitPM(n)
    cont <- 0
    maxoper <- 500
    for (i in 1:(n - 1)) {
        flag_ii_div = T
        while (flag_ii_div) {
            flag_ii_z = T
            while (flag_ii_z) {
                flag_row = T
                while (flag_row) {
                  nzcount = 0
                  for (j in (i + 1):n) {
                    if (abs_r(A[i, j]) > tol) {
                      nzcount = nzcount + 1
                      q = A[i, j]/A[i, i]
                      r = A[i, j]\%\%A[i, i]
                      if (abs_r(r) < tol) {
                        T1 = elemOperA(n, i, j, -q)
                        T1inv = elemOperA(n, i, j, q)
                        A = trimPolyMat(A, tol2)
                        A_ant = A
                        A = A * T1
                        U = T1inv * U
                        cont <- cont + 1
                        if (cont > maxoper) {
                          stop("smith_polymath:MaxOperExceeded", 
                            "Max Oper Exceeded")
                        }
                      }
                      else {
                        T1 = elemOperA(n, i, j, -q)
                        T1inv = elemOperA(n, i, j, q)
                        T2 = elemOperS(n, i, j)
                        A_ant = A
                        A = A * T1 * T2
                        A = trimPolyMat(A, tol2)
                        U = t(T2) * T1inv * U
                        cont = cont + 1
                        if (cont > maxoper) {
                          stop("smith_polymath:MaxOperExceeded", 
                            "Max Oper Exceeded")
                        }
                      }
                    }
                  }
                  if (nzcount == 0) {
                    flag_row = F
                  }
                }
                flag_col = T
                while (flag_col) {
                  nzcount = 0
                  for (j in (i + 1):n) {
                    if (abs_r(A[j, i]) > tol) {
                      nzcount = nzcount + 1
                      q = A[j, i]/A[i, i]
                      r = A[j, i]\%\%A[i, i]
                      if (abs_r(r) < tol) {
                        T1 = elemOperA(n, j, i, -q)
                        T1inv = elemOperA(n, j, i, q)
                        A_ant = A
                        A = T1 * A
                        A = trimPolyMat(A, tol2)
                        V = V * T1inv
                        cont = cont + 1
                        if (cont > maxoper) {
                          stop("smith_polymath:MaxOperExceeded", 
                            "Max Oper Exceeded")
                        }
                      }
                      else {
                        T1 = elemOperA(n, j, i, -q)
                        T1inv = elemOperA(n, j, i, q)
                        T2 = elemOperS(n, i, j)
                        A_ant = A
                        A = T2 * T1 * A
                        A = trimPolyMat(A, tol2)
                        V = V * T1inv * t(T2)
                        cont = cont + 1
                        if (cont > maxoper) {
                          stop("smith_polymath:MaxOperExceeded", 
                            "Max Oper Exceeded")
                        }
                      }
                    }
                  }
                  if (nzcount == 0) {
                    flag_col = F
                  }
                }
                flag_ii_z = F
                for (j in (i + 1):n) {
                  if (abs_r(A[i, j]) > tol | abs_r(A[j, i]) > 
                    tol) {
                    flag_ii_z = T
                  }
                }
            }
            flag_ii_div = F
            for (j in (i + 1):n) {
                for (k in (i + 1):n) {
                  q = A[j, k]/A[i, i]
                  r = A[j, k]\%\%A[i, i]
                  if (abs_r(r) > tol) {
                    flag_ii_div = T
                    T1 = elemOperA(n, i, j, polynomial(coef = c(1)))
                    T1inv = elemOperA(n, i, j, polynomial(coef = c(-11)))
                    A_ant = A
                    A = T1 * A
                    V = V * T1inv
                    cont = cont + 1
                    if (cont > maxoper) {
                      stop("smith_polymath:MaxOperExceeded", 
                        "Max Oper Exceeded")
                    }
                    break
                  }
                }
                if (flag_ii_div == 1) {
                  break
                }
            }
        }
    }
    D <- A
    return(list(U = U, D = D, V = V))
  }
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
\keyword{ ~kwd1 }
\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line
